########第一题#####
# 对数化所有蛋白数据
#第三次实习
pro1=read.table("ADdata1.txt",row.names = NULL)#读取文件，不将第一列作为行名
pro2=read.csv("ADdata2.csv")
library(openxlsx)
pro3=read.xlsx("ADdata3.xlsx")
pro4=read.table("ADdata4.txt",row.names = NULL)
h1=merge(pro1,pro2,by=1)#两两按第一列合并，最后再合并
h2=merge(pro3,pro4,by=1)
result=merge(h1,h2,by=1)
data_log2<-log2(result[,2:65])
# 提取分组信息
col_name<-colnames(data_log2)
asym<-grep("^asym",col_name)
ctl<-grep("^ctl",col_name)
ad<-grep("^ad",col_name)
data<-data_log2[,c(asym,ctl,ad)]
group<-c(rep("asym",15),rep("ctl",10),rep("ad",39))
result1<-data.frame(ID=result$row.names)
#选第二个蛋白尝试
pro<-unlist(data[93,])
num1<-sum(pro[asym]!=-Inf)
num2<-sum(pro[ctl]!=-Inf)
num3<-sum(pro[ad]!=-Inf)
 if(num1>2 & num2>2 & num3>2)
{
  test1<-oneway.test(pro~group)
  if(test1$p.value == "NaN")
    test1$p.value=NA
 }
print(test1$p.value)
# 循环-判断各组数据是否有至少3个样本是有数值的（非0或非NA）  
for(i in 1:100)
{
  pro<-unlist(data[i,])
  num1<-sum(pro[asym]!=-Inf)
  num2<-sum(pro[ctl]!=-Inf)
  num3<-sum(pro[ad]!=-Inf)
  if(num1 > 2 & num2 > 2 & num3 > 2)
  {
    test<-oneway.test(pro~group)
    if(test$p.value == "NaN")
      test$p.value=NA
  } 
  else
    test$p.value=NA
  result1$p_value[i]<- test$p.value
}

#######第二题######
#加载数据
load("volcano.RData") 
#筛选p < 0.05的蛋白
go_pro<- subset(prostat,P< 0.05)
#安装clusterProfiler包
if (!requireNamespace("BiocManager", quietly = TRUE))  
  install.packages("BiocManager")  
BiocManager::install("clusterProfiler")  
BiocManager::install("org.Hs.eg.db")
library(clusterProfiler)
library(org.Hs.eg.db)
#使用enrichGO函数进行GO富集分析
geneList <- go_pro$ID  # volcano_data中有一个列名为ID的列，存储了蛋白的ID  
ego <- enrichGO(gene         = geneList,   
                OrgDb        = 'org.Hs.eg.db', # 假设是人类的蛋白数据，选择合适的OrgDb  
                keyType      = 'SYMBOL',   
                universe     = NULL,   
                ont          = "BP",  # 可以是"BP" (生物过程), "CC" (细胞成分), 或 "MF" (分子功能)  
                pAdjustMethod = "none",   
                pvalueCutoff = 0.1,   
                qvalueCutoff = 0.2,   
                readable     = TRUE,
                minGSSize = 10,maxGSSize = 500,
                )
picture2<-dotplot(ego, showCategory=20, title="GO Enrichment Analysis")
print(picture2)
save(picture2,file="8305211202-何加星-实习五气泡图.jpg")
